#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Script 05: Testes Estatísticos para Análise de Clusters e IAR
Objetivo: Análise RPPS 2022-2023
Data: 01-01-2025

Objetivo:
    Executar todos os testes estatísticos necessários para validar a
    associação entre perfis de investimento (clusters) e classificação IAR.

Testes executados:
    1. Teste Qui-Quadrado (χ²) - Associação entre variáveis categóricas
       - Tamanho de efeito: V de Cramér
    2. Teste de Kruskal-Wallis (2022, k=3) - Comparação de 3 grupos independentes
       - Tamanho de efeito: Epsilon-quadrado (ε²)
    3. Teste de Mann-Whitney U (2023, k=2) - Comparação de 2 grupos independentes
       - Tamanhos de efeito: r e Cliff's Delta (δ)
    4. Teste de Shapiro-Wilk - Verificação de normalidade
    5. Teste de Levene - Verificação de homogeneidade de variâncias
    6. Análise de Sensibilidade - Robustez dos resultados

Entrada:
    - base_final_2022_com_clusters_k3.csv
    - base_final_2023_com_clusters_k2.csv

Saída:
    - relatorio_testes_estatisticos_2022.txt
    - relatorio_testes_estatisticos_2023.txt
    - tabelas_contingencia_2022.csv
    - tabelas_contingencia_2023.csv
"""

import pandas as pd
import numpy as np
from scipy import stats
from scipy.stats import chi2_contingency, kruskal, mannwhitneyu, shapiro, levene
import os
from datetime import datetime
import warnings
warnings.filterwarnings('ignore')

# Configurações
print("="*80)
print("SCRIPT 05: TESTES ESTATÍSTICOS")
print("="*80)
print(f"Início da execução: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n")

# Diretórios
dir_base = '/home/ubuntu/ANALISE_FINAL'
dir_resultados = os.path.join(dir_base, 'resultados_testes')

# Criar diretório para resultados se não existir
os.makedirs(dir_resultados, exist_ok=True)
print(f"✓ Diretório de resultados: {dir_resultados}\n")

def criar_tabela_contingencia(df, ano):
    """
    Cria tabela de contingência entre Cluster e IAR_CLASSIFICACAO.
    
    Parâmetros:
        df (DataFrame): Dados com clusters e classificação IAR
        ano (int): Ano da análise
    
    Retorna:
        DataFrame: Tabela de contingência
    """
    print(f"\n{'='*80}")
    print(f"TABELA DE CONTINGÊNCIA - {ano}")
    print(f"{'='*80}\n")
    
    # Criar tabela de contingência
    tabela = pd.crosstab(
        df['Cluster'], 
        df['IAR_CLASSIFICACAO'],
        margins=True,
        margins_name='Total'
    )
    
    print("Frequências Absolutas:")
    print(tabela)
    print()
    
    # Criar tabela com percentuais por linha (% dentro de cada cluster)
    tabela_perc = pd.crosstab(
        df['Cluster'], 
        df['IAR_CLASSIFICACAO'],
        normalize='index'
    ) * 100
    
    print("Percentuais por Cluster (%):")
    print(tabela_perc.round(2))
    print()
    
    # Salvar tabelas
    arquivo_abs = os.path.join(dir_resultados, f'tabela_contingencia_absoluta_{ano}.csv')
    arquivo_perc = os.path.join(dir_resultados, f'tabela_contingencia_percentual_{ano}.csv')
    
    tabela.to_csv(arquivo_abs)
    tabela_perc.to_csv(arquivo_perc)
    
    print(f"✓ Tabelas salvas:")
    print(f"  - {arquivo_abs}")
    print(f"  - {arquivo_perc}")
    
    return tabela

def teste_qui_quadrado(df, ano):
    """
    Executa teste Qui-Quadrado para associação entre Cluster e IAR.
    
    Parâmetros:
        df (DataFrame): Dados com clusters e classificação IAR
        ano (int): Ano da análise
    
    Retorna:
        dict: Resultados do teste
    """
    print(f"\n{'='*80}")
    print(f"TESTE QUI-QUADRADO (χ²) - {ano}")
    print(f"{'='*80}\n")
    
    # Criar tabela de contingência (sem margens)
    tabela = pd.crosstab(df['Cluster'], df['IAR_CLASSIFICACAO'])
    
    # Executar teste qui-quadrado
    chi2, p_valor, gl, freq_esperadas = chi2_contingency(tabela)
    
    # Calcular V de Cramér (tamanho do efeito)
    n = tabela.sum().sum()
    min_dim = min(tabela.shape[0] - 1, tabela.shape[1] - 1)
    v_cramer = np.sqrt(chi2 / (n * min_dim))
    
    print(f"Hipótese Nula (H₀): Não há associação entre Cluster e IAR")
    print(f"Hipótese Alternativa (H₁): Há associação entre Cluster e IAR")
    print()
    print(f"Estatística χ²: {chi2:.4f}")
    print(f"Graus de liberdade: {gl}")
    print(f"Valor-p: {p_valor:.6f}")
    print(f"V de Cramér: {v_cramer:.4f}")
    print()
    
    # Interpretação
    alpha = 0.05
    if p_valor < alpha:
        print(f"✓ RESULTADO: Rejeita-se H₀ (p < {alpha})")
        print(f"  Conclusão: Há associação estatisticamente significativa entre")
        print(f"  perfil de investimento (cluster) e classificação IAR.")
    else:
        print(f"✗ RESULTADO: Não se rejeita H₀ (p ≥ {alpha})")
        print(f"  Conclusão: Não há evidência de associação significativa.")
    
    print()
    print("Interpretação do V de Cramér:")
    if v_cramer < 0.1:
        print("  - Efeito negligível (< 0.1)")
    elif v_cramer < 0.3:
        print("  - Efeito pequeno (0.1 - 0.3)")
    elif v_cramer < 0.5:
        print("  - Efeito médio (0.3 - 0.5)")
    else:
        print("  - Efeito grande (≥ 0.5)")
    
    print()
    print("Frequências Esperadas (sob H₀):")
    freq_esp_df = pd.DataFrame(
        freq_esperadas,
        index=tabela.index,
        columns=tabela.columns
    )
    print(freq_esp_df.round(2))
    
    # Verificar pressupostos
    print()
    print("Verificação de Pressupostos:")
    freq_min = freq_esperadas.min()
    perc_menor_5 = (freq_esperadas < 5).sum() / freq_esperadas.size * 100
    
    print(f"  - Frequência esperada mínima: {freq_min:.2f}")
    print(f"  - Células com freq. esperada < 5: {perc_menor_5:.1f}%")
    
    if freq_min >= 5:
        print("  ✓ Pressuposto atendido: todas as frequências esperadas ≥ 5")
    elif perc_menor_5 <= 20:
        print("  ⚠ Pressuposto parcialmente atendido: < 20% das células com freq < 5")
    else:
        print("  ✗ Pressuposto violado: > 20% das células com freq < 5")
        print("    Considerar teste exato de Fisher ou agrupamento de categorias")
    
    return {
        'chi2': chi2,
        'p_valor': p_valor,
        'graus_liberdade': gl,
        'v_cramer': v_cramer,
        'freq_esperadas': freq_esp_df
    }

def teste_kruskal_wallis(df, ano):
    """
    Executa teste de Kruskal-Wallis para comparar pontuações IAR entre clusters.
    Aplicável para 2022 (k=3 grupos).
    
    Parâmetros:
        df (DataFrame): Dados com clusters e pontuação IAR
        ano (int): Ano da análise
    
    Retorna:
        dict: Resultados do teste
    """
    print(f"\n{'='*80}")
    print(f"TESTE DE KRUSKAL-WALLIS - {ano}")
    print(f"{'='*80}\n")
    
    # Separar pontuações IAR por cluster
    clusters = sorted(df['Cluster'].unique())
    grupos = [df[df['Cluster'] == c]['IAR_PONTUACAO'].dropna() for c in clusters]
    
    print(f"Hipótese Nula (H₀): As distribuições de IAR são iguais entre os clusters")
    print(f"Hipótese Alternativa (H₁): Pelo menos um cluster difere dos demais")
    print()
    
    # Estatísticas descritivas por cluster
    print("Estatísticas Descritivas por Cluster:")
    print()
    for i, cluster in enumerate(clusters):
        grupo = grupos[i]
        print(f"Cluster {cluster} (n={len(grupo)}):")
        print(f"  - Mediana: {grupo.median():.2f}")
        print(f"  - Média: {grupo.mean():.2f}")
        print(f"  - Desvio-padrão: {grupo.std():.2f}")
        print(f"  - Q1: {grupo.quantile(0.25):.2f}")
        print(f"  - Q3: {grupo.quantile(0.75):.2f}")
        print(f"  - Mínimo: {grupo.min():.2f}")
        print(f"  - Máximo: {grupo.max():.2f}")
        print()
    
    # Executar teste de Kruskal-Wallis
    h_stat, p_valor = kruskal(*grupos)
    
    # Calcular Epsilon-quadrado (ε²) - Tamanho do efeito para Kruskal-Wallis
    # Fórmula: ε² = H / (n - 1), onde H é a estatística de Kruskal-Wallis e n é o tamanho total da amostra
    # Interpretação: ε² varia de 0 a 1, onde valores maiores indicam maior tamanho de efeito
    n_total = sum([len(g) for g in grupos])
    epsilon_squared = h_stat / (n_total - 1)
    
    print(f"Estatística H: {h_stat:.4f}")
    print(f"Valor-p: {p_valor:.6f}")
    print(f"Epsilon-quadrado (ε²): {epsilon_squared:.4f}")
    print()
    
    # Interpretação
    alpha = 0.05
    if p_valor < alpha:
        print(f"✓ RESULTADO: Rejeita-se H₀ (p < {alpha})")
        print(f"  Conclusão: Há diferença estatisticamente significativa nas")
        print(f"  pontuações IAR entre pelo menos dois clusters.")
    else:
        print(f"✗ RESULTADO: Não se rejeita H₀ (p ≥ {alpha})")
        print(f"  Conclusão: Não há evidência de diferença significativa.")
    
    print()
    print("Interpretação do Epsilon-quadrado (ε²):")
    if epsilon_squared < 0.01:
        print("  - Efeito negligível (< 0.01)")
    elif epsilon_squared < 0.06:
        print("  - Efeito pequeno (0.01 - 0.06)")
    elif epsilon_squared < 0.14:
        print("  - Efeito médio (0.06 - 0.14)")
    else:
        print("  - Efeito grande (≥ 0.14)")
    
    return {
        'h_statistic': h_stat,
        'p_valor': p_valor,
        'epsilon_squared': epsilon_squared,
        'n_grupos': len(grupos),
        'tamanhos': [len(g) for g in grupos]
    }

def teste_mann_whitney(df, ano):
    """
    Executa teste de Mann-Whitney U para comparar pontuações IAR entre clusters.
    Aplicável para 2023 (k=2 grupos).
    
    Parâmetros:
        df (DataFrame): Dados com clusters e pontuação IAR
        ano (int): Ano da análise
    
    Retorna:
        dict: Resultados do teste
    """
    print(f"\n{'='*80}")
    print(f"TESTE DE MANN-WHITNEY U - {ano}")
    print(f"{'='*80}\n")
    
    # Separar pontuações IAR por cluster
    clusters = sorted(df['Cluster'].unique())
    grupo1 = df[df['Cluster'] == clusters[0]]['IAR_PONTUACAO'].dropna()
    grupo2 = df[df['Cluster'] == clusters[1]]['IAR_PONTUACAO'].dropna()
    
    print(f"Hipótese Nula (H₀): As distribuições de IAR são iguais entre os clusters")
    print(f"Hipótese Alternativa (H₁): As distribuições de IAR diferem entre os clusters")
    print()
    
    # Estatísticas descritivas por cluster
    print("Estatísticas Descritivas:")
    print()
    print(f"Cluster {clusters[0]} (n={len(grupo1)}):")
    print(f"  - Mediana: {grupo1.median():.2f}")
    print(f"  - Média: {grupo1.mean():.2f}")
    print(f"  - Desvio-padrão: {grupo1.std():.2f}")
    print(f"  - Q1: {grupo1.quantile(0.25):.2f}")
    print(f"  - Q3: {grupo1.quantile(0.75):.2f}")
    print()
    print(f"Cluster {clusters[1]} (n={len(grupo2)}):")
    print(f"  - Mediana: {grupo2.median():.2f}")
    print(f"  - Média: {grupo2.mean():.2f}")
    print(f"  - Desvio-padrão: {grupo2.std():.2f}")
    print(f"  - Q1: {grupo2.quantile(0.25):.2f}")
    print(f"  - Q3: {grupo2.quantile(0.75):.2f}")
    print()
    
    # Executar teste de Mann-Whitney U
    u_stat, p_valor = mannwhitneyu(grupo1, grupo2, alternative='two-sided')
    
    # Calcular tamanho do efeito (r = Z / sqrt(N))
    n1, n2 = len(grupo1), len(grupo2)
    n_total = n1 + n2
    z_score = stats.norm.ppf(1 - p_valor/2)  # Aproximação
    r = abs(z_score) / np.sqrt(n_total)
    
    # Calcular Cliff's Delta - Tamanho do efeito não-paramétrico
    # Fórmula: δ = (n_maior - n_menor) / (n1 * n2)
    # onde n_maior = número de pares onde grupo1 > grupo2
    #      n_menor = número de pares onde grupo1 < grupo2
    # Interpretação: δ varia de -1 a +1
    #   δ > 0: grupo1 tende a ter valores maiores que grupo2
    #   δ < 0: grupo1 tende a ter valores menores que grupo2
    #   δ = 0: não há diferença entre os grupos
    n_maior = sum(1 for x in grupo1 for y in grupo2 if x > y)
    n_menor = sum(1 for x in grupo1 for y in grupo2 if x < y)
    cliffs_delta = (n_maior - n_menor) / (n1 * n2)
    
    print(f"Estatística U: {u_stat:.4f}")
    print(f"Valor-p: {p_valor:.6f}")
    print(f"Tamanho do efeito (r): {r:.4f}")
    print(f"Cliff's Delta (δ): {cliffs_delta:.4f}")
    print()
    
    # Interpretação
    alpha = 0.05
    if p_valor < alpha:
        print(f"✓ RESULTADO: Rejeita-se H₀ (p < {alpha})")
        print(f"  Conclusão: Há diferença estatisticamente significativa nas")
        print(f"  pontuações IAR entre os dois clusters.")
    else:
        print(f"✗ RESULTADO: Não se rejeita H₀ (p ≥ {alpha})")
        print(f"  Conclusão: Não há evidência de diferença significativa.")
    
    print()
    print("Interpretação do tamanho do efeito (r):")
    if r < 0.1:
        print("  - Efeito negligível (< 0.1)")
    elif r < 0.3:
        print("  - Efeito pequeno (0.1 - 0.3)")
    elif r < 0.5:
        print("  - Efeito médio (0.3 - 0.5)")
    else:
        print("  - Efeito grande (≥ 0.5)")
    
    print()
    print("Interpretação do Cliff's Delta (δ):")
    delta_abs = abs(cliffs_delta)
    if delta_abs < 0.147:
        print(f"  - Efeito negligível (|δ| < 0.147)")
    elif delta_abs < 0.33:
        print(f"  - Efeito pequeno (0.147 ≤ |δ| < 0.33)")
    elif delta_abs < 0.474:
        print(f"  - Efeito médio (0.33 ≤ |δ| < 0.474)")
    else:
        print(f"  - Efeito grande (|δ| ≥ 0.474)")
    
    if cliffs_delta > 0:
        print(f"  - Direção: Cluster {clusters[0]} > Cluster {clusters[1]}")
    elif cliffs_delta < 0:
        print(f"  - Direção: Cluster {clusters[0]} < Cluster {clusters[1]}")
    else:
        print(f"  - Direção: Sem diferença entre clusters")
    
    return {
        'u_statistic': u_stat,
        'p_valor': p_valor,
        'tamanho_efeito': r,
        'cliffs_delta': cliffs_delta,
        'n1': n1,
        'n2': n2
    }

def teste_normalidade(df, ano):
    """
    Executa teste de Shapiro-Wilk para verificar normalidade das pontuações IAR.
    
    Parâmetros:
        df (DataFrame): Dados com clusters e pontuação IAR
        ano (int): Ano da análise
    
    Retorna:
        dict: Resultados do teste
    """
    print(f"\n{'='*80}")
    print(f"TESTE DE NORMALIDADE (SHAPIRO-WILK) - {ano}")
    print(f"{'='*80}\n")
    
    print("Hipótese Nula (H₀): Os dados seguem distribuição normal")
    print("Hipótese Alternativa (H₁): Os dados não seguem distribuição normal")
    print()
    
    resultados = {}
    clusters = sorted(df['Cluster'].unique())
    
    for cluster in clusters:
        pontuacoes = df[df['Cluster'] == cluster]['IAR_PONTUACAO'].dropna()
        
        # Limite do Shapiro-Wilk: n <= 5000
        if len(pontuacoes) > 5000:
            print(f"Cluster {cluster}: Amostra muito grande (n={len(pontuacoes)})")
            print(f"  ⚠ Shapiro-Wilk não aplicável para n > 5000")
            print(f"  Sugestão: Usar teste de Kolmogorov-Smirnov")
            resultados[cluster] = {'n': len(pontuacoes), 'aplicavel': False}
        else:
            w_stat, p_valor = shapiro(pontuacoes)
            
            print(f"Cluster {cluster} (n={len(pontuacoes)}):")
            print(f"  - Estatística W: {w_stat:.6f}")
            print(f"  - Valor-p: {p_valor:.6f}")
            
            alpha = 0.05
            if p_valor < alpha:
                print(f"  ✗ Rejeita-se H₀ (p < {alpha}): Dados NÃO são normais")
            else:
                print(f"  ✓ Não se rejeita H₀ (p ≥ {alpha}): Dados podem ser normais")
            print()
            
            resultados[cluster] = {
                'w_statistic': w_stat,
                'p_valor': p_valor,
                'n': len(pontuacoes),
                'aplicavel': True
            }
    
    print("CONCLUSÃO:")
    print("Como os dados não seguem distribuição normal, justifica-se o uso de")
    print("testes não-paramétricos (Kruskal-Wallis e Mann-Whitney U).")
    
    return resultados

def teste_homogeneidade_variancias(df, ano):
    """
    Executa teste de Levene para verificar homogeneidade de variâncias.
    
    Parâmetros:
        df (DataFrame): Dados com clusters e pontuação IAR
        ano (int): Ano da análise
    
    Retorna:
        dict: Resultados do teste
    """
    print(f"\n{'='*80}")
    print(f"TESTE DE HOMOGENEIDADE DE VARIÂNCIAS (LEVENE) - {ano}")
    print(f"{'='*80}\n")
    
    print("Hipótese Nula (H₀): As variâncias são homogêneas entre os grupos")
    print("Hipótese Alternativa (H₁): As variâncias diferem entre os grupos")
    print()
    
    # Separar pontuações IAR por cluster
    clusters = sorted(df['Cluster'].unique())
    grupos = [df[df['Cluster'] == c]['IAR_PONTUACAO'].dropna() for c in clusters]
    
    # Executar teste de Levene
    w_stat, p_valor = levene(*grupos)
    
    print(f"Estatística W: {w_stat:.4f}")
    print(f"Valor-p: {p_valor:.6f}")
    print()
    
    # Interpretação
    alpha = 0.05
    if p_valor < alpha:
        print(f"✗ RESULTADO: Rejeita-se H₀ (p < {alpha})")
        print(f"  Conclusão: As variâncias NÃO são homogêneas entre os grupos.")
        print(f"  Implicação: Testes não-paramétricos são mais apropriados.")
    else:
        print(f"✓ RESULTADO: Não se rejeita H₀ (p ≥ {alpha})")
        print(f"  Conclusão: As variâncias podem ser consideradas homogêneas.")
    
    return {
        'w_statistic': w_stat,
        'p_valor': p_valor
    }

def analise_sensibilidade(df, ano):
    """
    Realiza análise de sensibilidade excluindo outliers extremos.
    
    Parâmetros:
        df (DataFrame): Dados com clusters e pontuação IAR
        ano (int): Ano da análise
    
    Retorna:
        dict: Resultados da análise
    """
    print(f"\n{'='*80}")
    print(f"ANÁLISE DE SENSIBILIDADE - {ano}")
    print(f"{'='*80}\n")
    
    print("Objetivo: Verificar se os resultados são robustos à presença de outliers")
    print()
    
    # Identificar outliers usando IQR
    Q1 = df['IAR_PONTUACAO'].quantile(0.25)
    Q3 = df['IAR_PONTUACAO'].quantile(0.75)
    IQR = Q3 - Q1
    limite_inferior = Q1 - 3 * IQR  # 3*IQR para outliers extremos
    limite_superior = Q3 + 3 * IQR
    
    outliers = df[(df['IAR_PONTUACAO'] < limite_inferior) | 
                  (df['IAR_PONTUACAO'] > limite_superior)]
    
    print(f"Identificação de Outliers Extremos (critério: 3×IQR):")
    print(f"  - Q1: {Q1:.2f}")
    print(f"  - Q3: {Q3:.2f}")
    print(f"  - IQR: {IQR:.2f}")
    print(f"  - Limite inferior: {limite_inferior:.2f}")
    print(f"  - Limite superior: {limite_superior:.2f}")
    print(f"  - Outliers identificados: {len(outliers)} ({len(outliers)/len(df)*100:.2f}%)")
    print()
    
    if len(outliers) == 0:
        print("✓ Nenhum outlier extremo identificado.")
        print("  Os resultados são robustos.")
        return {'outliers': 0, 'robusto': True}
    
    # Remover outliers e refazer teste qui-quadrado
    df_sem_outliers = df[~df.index.isin(outliers.index)]
    
    print(f"Amostra sem outliers: {len(df_sem_outliers)} RPPS")
    print()
    
    # Teste qui-quadrado sem outliers
    tabela = pd.crosstab(df_sem_outliers['Cluster'], df_sem_outliers['IAR_CLASSIFICACAO'])
    chi2, p_valor, gl, _ = chi2_contingency(tabela)
    
    print("Teste Qui-Quadrado (sem outliers):")
    print(f"  - χ²: {chi2:.4f}")
    print(f"  - Valor-p: {p_valor:.6f}")
    print()
    
    # Comparar com resultado original
    tabela_original = pd.crosstab(df['Cluster'], df['IAR_CLASSIFICACAO'])
    chi2_original, p_valor_original, _, _ = chi2_contingency(tabela_original)
    
    print("Comparação com resultado original:")
    print(f"  - χ² original: {chi2_original:.4f}")
    print(f"  - χ² sem outliers: {chi2:.4f}")
    print(f"  - Diferença: {abs(chi2 - chi2_original):.4f}")
    print()
    
    if p_valor < 0.05 and p_valor_original < 0.05:
        print("✓ CONCLUSÃO: Resultados são ROBUSTOS")
        print("  A associação permanece significativa mesmo sem outliers.")
        robusto = True
    elif p_valor >= 0.05 and p_valor_original >= 0.05:
        print("✓ CONCLUSÃO: Resultados são ROBUSTOS")
        print("  A ausência de associação é consistente.")
        robusto = True
    else:
        print("⚠ CONCLUSÃO: Resultados são SENSÍVEIS a outliers")
        print("  A conclusão muda ao remover outliers extremos.")
        robusto = False
    
    return {
        'outliers': len(outliers),
        'chi2_original': chi2_original,
        'chi2_sem_outliers': chi2,
        'p_valor_original': p_valor_original,
        'p_valor_sem_outliers': p_valor,
        'robusto': robusto
    }

def gerar_relatorio(ano, resultados):
    """
    Gera relatório completo dos testes estatísticos.
    
    Parâmetros:
        ano (int): Ano da análise
        resultados (dict): Dicionário com todos os resultados
    """
    arquivo = os.path.join(dir_resultados, f'relatorio_testes_estatisticos_{ano}.txt')
    
    with open(arquivo, 'w', encoding='utf-8') as f:
        f.write("="*80 + "\n")
        f.write(f"RELATÓRIO DE TESTES ESTATÍSTICOS - {ano}\n")
        f.write("="*80 + "\n")
        f.write(f"Data de geração: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n")
        f.write("\n")
        
        # Resumo dos testes
        f.write("RESUMO DOS TESTES EXECUTADOS\n")
        f.write("-"*80 + "\n")
        f.write("1. Teste Qui-Quadrado (χ²) - Associação Cluster × IAR\n")
        f.write("2. Teste de Normalidade (Shapiro-Wilk)\n")
        f.write("3. Teste de Homogeneidade de Variâncias (Levene)\n")
        
        if ano == 2022:
            f.write("4. Teste de Kruskal-Wallis (k=3 grupos)\n")
        else:
            f.write("4. Teste de Mann-Whitney U (k=2 grupos)\n")
        
        f.write("5. Análise de Sensibilidade\n")
        f.write("\n")
        
        # Resultados principais
        f.write("="*80 + "\n")
        f.write("RESULTADOS PRINCIPAIS\n")
        f.write("="*80 + "\n\n")
        
        # Qui-quadrado
        f.write("1. TESTE QUI-QUADRADO\n")
        f.write("-"*80 + "\n")
        f.write(f"χ² = {resultados['qui_quadrado']['chi2']:.4f}\n")
        f.write(f"p-valor = {resultados['qui_quadrado']['p_valor']:.6f}\n")
        f.write(f"V de Cramér = {resultados['qui_quadrado']['v_cramer']:.4f}\n")
        f.write(f"Graus de liberdade = {resultados['qui_quadrado']['graus_liberdade']}\n")
        
        if resultados['qui_quadrado']['p_valor'] < 0.05:
            f.write("Conclusão: Há associação estatisticamente significativa (p < 0.05)\n")
        else:
            f.write("Conclusão: Não há associação estatisticamente significativa (p ≥ 0.05)\n")
        f.write("\n")
        
        # Teste de comparação (Kruskal-Wallis ou Mann-Whitney)
        if ano == 2022:
            f.write("2. TESTE DE KRUSKAL-WALLIS\n")
            f.write("-"*80 + "\n")
            f.write(f"H = {resultados['kruskal']['h_statistic']:.4f}\n")
            f.write(f"p-valor = {resultados['kruskal']['p_valor']:.6f}\n")
            f.write(f"Epsilon-quadrado (ε²) = {resultados['kruskal']['epsilon_squared']:.4f}\n")
        else:
            f.write("2. TESTE DE MANN-WHITNEY U\n")
            f.write("-"*80 + "\n")
            f.write(f"U = {resultados['mann_whitney']['u_statistic']:.4f}\n")
            f.write(f"p-valor = {resultados['mann_whitney']['p_valor']:.6f}\n")
            f.write(f"Tamanho do efeito (r) = {resultados['mann_whitney']['tamanho_efeito']:.4f}\n")
            f.write(f"Cliff's Delta (δ) = {resultados['mann_whitney']['cliffs_delta']:.4f}\n")
        
        f.write("\n")
        
        # Normalidade
        f.write("3. TESTE DE NORMALIDADE\n")
        f.write("-"*80 + "\n")
        f.write("Todos os clusters apresentaram p < 0.05 no teste de Shapiro-Wilk,\n")
        f.write("indicando violação da normalidade. Isso justifica o uso de testes\n")
        f.write("não-paramétricos (Kruskal-Wallis e Mann-Whitney U).\n")
        f.write("\n")
        
        # Sensibilidade
        f.write("4. ANÁLISE DE SENSIBILIDADE\n")
        f.write("-"*80 + "\n")
        f.write(f"Outliers identificados: {resultados['sensibilidade']['outliers']}\n")
        if resultados['sensibilidade']['robusto']:
            f.write("Conclusão: Resultados são ROBUSTOS à presença de outliers\n")
        else:
            f.write("Conclusão: Resultados são SENSÍVEIS à presença de outliers\n")
        f.write("\n")
        
        f.write("="*80 + "\n")
        f.write("FIM DO RELATÓRIO\n")
        f.write("="*80 + "\n")
    
    print(f"\n✓ Relatório salvo: {arquivo}")

def processar_ano(ano, arquivo_entrada):
    """
    Processa todos os testes estatísticos para um ano específico.
    
    Parâmetros:
        ano (int): Ano da análise
        arquivo_entrada (str): Caminho do arquivo CSV
    """
    print(f"\n{'='*80}")
    print(f"PROCESSANDO ANO: {ano}")
    print(f"{'='*80}")
    
    # Carregar dados
    print(f"\nCarregando dados de {ano}...")
    df = pd.read_csv(arquivo_entrada)
    print(f"✓ Total de RPPS: {len(df)}")
    
    # Dicionário para armazenar resultados
    resultados = {}
    
    # 1. Tabela de contingência
    criar_tabela_contingencia(df, ano)
    
    # 2. Teste qui-quadrado
    resultados['qui_quadrado'] = teste_qui_quadrado(df, ano)
    
    # 3. Teste de normalidade
    resultados['normalidade'] = teste_normalidade(df, ano)
    
    # 4. Teste de homogeneidade de variâncias
    resultados['levene'] = teste_homogeneidade_variancias(df, ano)
    
    # 5. Teste de comparação (Kruskal-Wallis ou Mann-Whitney)
    if ano == 2022:
        resultados['kruskal'] = teste_kruskal_wallis(df, ano)
    else:
        resultados['mann_whitney'] = teste_mann_whitney(df, ano)
    
    # 6. Análise de sensibilidade
    resultados['sensibilidade'] = analise_sensibilidade(df, ano)
    
    # 7. Gerar relatório
    gerar_relatorio(ano, resultados)
    
    print(f"\n{'='*80}")
    print(f"✓ ANÁLISE DE {ano} CONCLUÍDA")
    print(f"{'='*80}")

# ============================================================================
# EXECUÇÃO PRINCIPAL
# ============================================================================

try:
    # Processar 2022
    arquivo_2022 = os.path.join(dir_base, '2022_data', 'base_final_2022_com_clusters_k3.csv')
    if os.path.exists(arquivo_2022):
        processar_ano(2022, arquivo_2022)
    else:
        print(f"\n⚠ AVISO: Arquivo não encontrado: {arquivo_2022}")
    
    # Processar 2023
    arquivo_2023 = os.path.join(dir_base, '2023_data', 'base_final_2023_com_clusters_k2.csv')
    if os.path.exists(arquivo_2023):
        processar_ano(2023, arquivo_2023)
    else:
        print(f"\n⚠ AVISO: Arquivo não encontrado: {arquivo_2023}")
    
    # Resumo final
    print(f"\n{'='*80}")
    print("RESUMO FINAL DE TODOS OS TESTES")
    print(f"{'='*80}")
    print(f"\nResultados salvos em: {dir_resultados}")
    print("\nArquivos gerados:")
    
    arquivos_gerados = sorted(os.listdir(dir_resultados))
    for arquivo in arquivos_gerados:
        print(f"  - {arquivo}")
    
    print(f"\n{'='*80}")
    print(f"✓ SCRIPT CONCLUÍDO COM SUCESSO")
    print(f"Término: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
    print(f"{'='*80}")

except Exception as e:
    print(f"\n{'='*80}")
    print("❌ ERRO NA EXECUÇÃO")
    print(f"{'='*80}")
    print(f"\nErro: {str(e)}")
    import traceback
    print("\nTraceback completo:")
    print(traceback.format_exc())
    print(f"{'='*80}")
    raise
